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A simple yet robust and accurate approach for capturing shock waves using a high-order 
discontinuous Galerkin (DG) method is presented. The method uses the physical viscous 
terms of the Navier-Stokes equations as suggested by others; however, the proposed for- 
mulation of the numerical viscosity is continuous and compact by construction, and does 
not require the solution of an auxiliary diffusion equation. This work also presents two 
analyses that guided the formulation of the numerical viscosity and certain aspects of the 
DG implementation. A local eigenvalue analysis of the DG discretization applied to a shock 
containing element is used to evaluate the robustness of several Riemann flux functions, 
and to evaluate algorithm choices that exist within the underlying DG discretization. A 
second analysis examines exact solutions to the DG discretization in a shock containing 
element, and identifies a “model” instability that will inevitably arise when solving the 
Euler equations using the DG method. This analysis identifies the minimum viscosity re- 
quired for stability. The shock capturing method is demonstrated for high-speed flow over 
an inviscid cylinder and for an unsteady disturbance in a hypersonic boundary layer. Nu- 
merical tests are presented that evaluate several aspects of the shock detection terms. The 
sensitivity of the results to model parameters is examined with grid and order refinement 
studies. 


I. Introduction 

The discontinuous Galerkin (DG) method is quickly becoming a popular means to obtain high-order 
solutions for a broad range of flows. However, flows with strong shocks continue to be a challenge for all high- 
order methods. Techniques for capturing shocks using DG methods range from simply reducing the order of 
the method 1-2 or adding viscosity in the vicinity of a shock, 3-6 to carefully constructed limiters. 7-11 Local 
order reduction, combined with aggressive mesh adaptation, may give acceptable results in many problems. 
However, it introduces many new complications into the software implementation: shock detection, mesh 
movement, mesh refinement, load rebalancing, and added computational expense, to mention a few. Limiter 
techniques are commonly used and have been demonstrated for moderate order (p=3). However, their 
formulations tend not to be generalized (e.g. tied to particular quadrature sets and element shapes), and 
they tend to become non-compact as the order increases. There are no clear or straight forward extensions 
to arbitrary element shapes and arbitrarily high-order for many limiter formulations. 

The method of artificially increasing the physical viscosity in the vicinity of a shock is attractive in that it 
is readily applied to any existing Navier-Stokes code with no reduction of the formal order. Shock capturing 
based on an artificial viscosity approach is commonly used with low order finite-difference and finite-volume 
methods. For many formulations, there is a direct link between flux limiter and artificial viscosity approaches. 
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Persson and Peraire 5 were the first to apply the technique to DG for solutions of polynomial degree much 
greater than 1. Their approach was to use sufficiently high order and sufficiently high viscosity such that 
the shock is resolved within the span of a single element. They formulated a piecewise constant artificial 
viscosity that depended directly on a local shock sensor. Later work by Barter and Darmofal 6 demonstrated 
the benefits of having a smooth numerical viscosity, which they formulated by solving an auxiliary diffusion 
equation along with the usual governing equations. Both works used strong implicit solvers which tend to 
suppress robustness issues that will frequently appear in large 3D problems when the same implicit solver is 
less effective or even impractical. The approach has been qualitatively demonstrated on a range of problems 
from ID shock-tubes, to transonic airfoils and high-Mach blunt body problems. 

This work presents a simple algebraic algorithm for determining the artificial viscosity in the vicinity 
of a shock. The approach produces a continuous field without the need to solve an additional partial 
differential equation. As such, it is simple to implement and adds little computational expense. The method 
is demonstrated for steady and unsteady shock waves in both viscous and inviscicl simulations. Away from 
a shock wave, where the flow is smooth, the viscosity decays at a rate of p + 2. 

Most efforts for the design and analysis of higlr-order methods focus on efficiency and asymptotic proper- 
ties. However, the viscous shock capturing approach requires that a captured shock will remain marginally 
resolved as the mesh is refined. Two types of analysis are presented that provide insight into the behavior of 
the DG discretization in the vicinity of a captured shock wave. This insight allows the DG implementation 
to be tailored to capture shock waves with minimal artificial viscosity. 

The first section of this article describes the Navier-Stokes equations and the application of DG to those 
equations. The next section describes two analyses and some observations drawn from them. A local 
eigenvalue analysis of inviscid flow in the vicinity of a shock containing element is used to evaluate several 
approximate Riemann flux methods, and to evaluate the effect of approximating volume and edge fluxes by 
various means. The second analysis presents exact solutions to the DG discretization of the Euler equations. 
This analysis reveals an inherent instability that must be mitigated through dissipation or other means. The 
third section introduces a compact procedure for constructing a continuous artificial viscosity that does not 
require solving an auxiliary partial-differntial equation. The last section presents numerical experiments that 
evaluate the performance and sensitivity of the method to user-supplied input and other factors. 


II. Governing Equations and DG Formulation 


A. Navier-Stokes Equations 

The non-dimensional compressible Navier-Stokes equations for a perfect gas in conservation form, using 
tensor index notation are: 


dp d{pu 3 ) = 

dt dxj 

dpe d ( huj - UjTjj + qj) _ 

dt dxj ’ 

dpui d (puiUj + SijP - Tij) 

dt dxj 


( 1 ) 

( 2 ) 

i = 1,2,3 (3) 


where p is the density, e is the internal energy per unit mass, P is the pressure, rq is the component of 
velocity in the Cartesian coordinate direction x,;, and 5ij is the Kronecker delta. The quantities h, qt, and 
Tij are the enthalpy, heat flux, and shear stress terms, respectively; and are given by: 


h = pe + P, 


-7 1 p dT 

^ 7 — 1 Pr Re , dxi ’ 

_ p, fduj d'Uj 2 du k \ 

Re r \dxj dxj l ’ j 3dxk) 

where Pr is the Prandtl number, and Re r is the Reynolds number based on reference state of the non- 
dimensionalization, and T is the tempterature given by T = P/ p = (7 — l)(e — itfcWfc/2). The length and 


2 of 22 


American Institute of Aeronautics and Astronautics 



the thermodynamic variables have been non-dimensionalized with respect to a prescribed reference state: 
Xi = Xi/L r , p = p/p r , P = P/P r , T = T/T r , and p = p/p r where denotes dimensional quantities, 
and the subscript r denotes the reference state. Other variables are normalized with respect to derived 

reference states as follows: u = u/u r , t = t/(L r /u r ), e = e r /u 2 , and Re r = p r u r L r /p r where u r = \Jp r /p r . 

Sutherland’s formula is used to evaluate both p r as a function of T r . and p, as a function of T. 

DG is applied to each of Eqs. 1-3 in essentially the same manner; however, the inviscid and viscous terms 
are treated differently. To facilitate the following discussion, each equation is cast in the general form: 

+ V • (Fj — F„) = 0, (4) 

where U denotes either p, pe, or piq, and Fj and F,, denote the inviscid and viscous contributions to the flux. 
Eq. (4) is transformed to a local computational coordinate system for each discrete element. Let (£, 77 , () 
denote the local coordinate system with Jacobian J = d(xi,X 2 ,X 3 )/d(£,r)X)- Equation (4) has the same 
form in the local coordinates: 

+ V • (Fj — F. u ) = 0, (5) 

but with U = U\J\, Fj = J -1 Fj|J|, and F v = J _1 F„| J|. 

The viscous shock capturing method involves augmenting the physical viscosity with a numerical artificial 
viscosity using, for example, p = p n + p P , where p„ denotes the numerical artificial viscosity, and p p denotes 
the physical viscosity. A compact method for computing the numerical viscosity is presented in section IV. 
One objective of the viscous shock capturing approach is to activate the numerical viscosity only in the 
vicinity of a shock wave. However, due to the heuristic nature of typical shock sensors, it is not unusual 
for the numerical viscosity to be non-zero in regions away from shock waves. This is especially a problem 
on coarse grids in which smooth flow features may be only marginally resolved, and may be mistaken for 
a shock wave. These difficulties can be minimized by using alternate means of introducing the artificial 
viscosity such as p = ma x(p n ,p p ), or by applying a flow dependent masking function that explicitly zeros 
the artificial viscosity in regions well away from where shock waves are know to exist. 


B. DG Discretization 

The implementation of DG used in this study follows the quadrature-free form described in Refs. 12-14. 
The following description primarily serves to introduce notation and identify terms that are the subject of 
the analyses and modifications related to the viscous shock capturing method described in later sections. 
Further details concerning the DG method in its implementation can be found in the references given above. 

The computational domain is subdivided into non-overlapping elements that cover the domain. The DG 
discretization is formulated locally in each element in a similar manner. The solution within each element is 
approximated as an expansion in a local basis set {b k }, usually polynomials of degree < p , 

N p 

U = ^2 b k u k , 

k = 0 

where N p denotes the number of terms in the basis set of degree p. The current work uses two-dimensional 
monomials of the form for all (i. j) pairs satisfying 0 < i + j < p. The number of unknowns in each 
element is the number of physical variables times the size of the basis set, N p . An equal number of equations 
governing these unknowns is derived by multiplying the governing equations by each member of the basis 
set, and integrating over the element. The integrals of the flux terms are integrated by parts to obtain the 
following weak form: 

[ b k ^dCl- [ V6fc(Fj — F„) dCl + [ 6 fe (Fj - F„) • n ds = 0 Mk : 0 < k < N p , (6) 

J n bit, J n J dn 

where fi denotes the element, dfl denotes the element boundary, and n denotes the outward unit normal. 
Following the quadrature-free formulation, 12 the fluxes are expanded in a similar manner, 

M M 

Fj = ^2 b k fi,k and F,. = ^ b k f v ,k- 
k — 0 k — 0 
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The degree of the flux expansions is allowed to be higher than that of the solution with the result that 
M > N p . The edge integral is evaluated piecewise on each discrete edge segment, <9%, defining the boundary 
of an individual element. Neighboring elements on either side of an edge segment have independent local 
approximations for the solution and fluxes. To resolve this ambiguity, the fluxes in the edge integrals are 
replaced by numerical fluxes that are a function of the solutions in both neighboring elements. The numerical 
edge fluxes are represented in terms of a lower dimensional edge basis, b k ,f 

M M 

F,; * n| af2 ^ s Fi — 'y ' bkjfi t k,j and F^ • n| 0fi ^ s F v = ^ ' b k jf v , k ,ji 
k = 0 k—0 

where dflj denotes an individual edge and M denotes the number of terms in the edge basis. 

The inviscid edge flux is modeled by an approximate Riemann flux. The local Lax-Friedrichs flux is 
simple and inexpensive to implement, and has worked well with DG for smooth flows. 14-1 ' As will be shown, 
however, it is not suitable for flow near shock waves. This work will analyze and compare the Roe flux 18 
and the Harten-Lax-van Leer (HLL) flux. 19 The Roe flux is well known for its robust ability to capture 
shocks when used in 2nd-order finite-volume methods. The HLL flux can also capture shocks well, and like 
the local Lax-Friedrichs flux, is very simple and efficient to implement. It is not popular among low-order 
finite-volume methods because of its inability to resolve contact discontinuities; however, several extensions 
to the basic HLL flux have been proposed 19,20 to remedy this weakness. Away from shock waves, the HLL 
flux performs similar to the Lax-Friedrichs flux and is sufficient for the present work. 

The solution gradients required for the viscous terms are evaluated using the DG methodology in a similar 
manner. Let 

N p 

<r w = ^2 b k crk = Vru, 

k = 0 

where w denotes either m or T. It immediately follows that 


b k a w dVL — / X7b k w dtt + / b k w n ds = 0. 


Id n 


As before, the multi-valued edge flux is replaced by a numerical flux, u>n| a£J = wn. Both F v and win are 
evaluated as described in Ref. 14. 

Because all of the fluxes are represented as expansions in the basis set, the integrations can be performed 
directly to produce a matrix equation of the form: 


M 


~du k 

dt 


+ V • [f^fc — + ^2 


fy fc,7 f v,k,j 


= o, 


where M, V and B are derived in Ref. 12 along with further details of the quadrature-free formulation. 


C. Non-linear Flux Terms 

The quadrature-free form of DG requires that all fluxes be expanded in the basis set. This process is trivial 
and exact when the equations are linear, but must be approximated for most non-linear cases. The inviscid 
fluxes of the Euler equations consist of ratios of polynomials such as 

jr = (. PUiKpuj ) Qr j = (pui){pe) 

P P 

As previously described in Refs. 12 and 13, polynomial approximations to the flux can be formulated by 
multiplying the fluxes by /?, and projecting back into the polynomial space. For example: 

J b k pf = j b k (pui)(puj) Vfc : 0 < k < M. (7) 

Considering only design order arguments, the polynomial products can be truncated to degree p, and the 
projection can be limited to include all k < M = N p -i. However, it was shown in Ref. 13 that retaining 
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the higher-order terms of polynomial products improved the robustness when capturing shocks. Polynomial 
products are not truncated in the current work, and the projection is onto the same polynomial space as the 
solution, k < M = N p . The operation is efficiently performed by noting that the left-hand side of Eq. 7 is 
the same for all terms, allowing the L 17-decomposition of the matrix to be re-used many times. 

The u t and T terms, needed to compute the viscous gradients, are obtained in a similar manner by 
solving: 

J bkpf = R V k : 0 < k < M, where R = j bkP or J bk{pui) i = 1,2, 3, 

and noting that the pressure, P, is a quadratic function of the dependent variables. 

III. Analysis of DG for a Shock Containing Element 

Two types of analysis are used to gain insight into the behavior of DG in the vicinity of a shock. The 
first analysis involves examining the eigenvalues of a DG discretization that has been linearized about a 
trial solution representing a shock. This analysis is an extension of an earlier study 13 that examined the 
DG method for the non-linear Burgers equation. There it was shown that simply retaining the higher-order 
terms of polynomial products resulted in a stable method for p < 3, without the need for limiters or added 
dissipation. For p > 3, the DG discretization with fully expanded products was unstable, but only for a 
small region of the solution space where the shock is very close to the edge of an element. This instability 
could be eliminated for Burgers equation with p < 5 by using a specially formulated approximate Riemann 
flux. The analysis is used here to examine the impact of different Riemann fluxes, of truncating polynomial 
products, and other aspects of the DG implementation. 

The second analysis takes advantage of special features of DG to construct exact steady solutions to 
the DG discretization in a shock containing element. Examining these solutions provides insight into the 
behavior of DG in such element. In particular, the analysis shows that the pressure can become negative 
leading to an instability. 

A. Eigenvalue analysis 

The eigenvalue analysis examines the DG discretization of the ID Euler equations for several elements in the 
vicinity of a shock. Given a trial solution representing a shock at a specified location within a unit element 
— | < x <§, the DG discretization is evaluated in the shock containing element and in two unit elements on 
each side. The trial solution is obtained by projecting the exact piecewise constant solution onto the discrete 
polynomial solution space. The resulting system of discrete equations is linearized about the trial solution 
and its eigenvalues are examined. The system is cast such that eigenvalues in the left half plane are stable. 
The analysis is carried out using the symbolic software Maple. 

The eigenvalue analysis is used to examine several choices for the approximate Riemann flux. Figure 1 
gives the largest real component of the eigenvalues as a function of shock position for a trial solution with 
a shock Mach number of 1.3, and for degree p ranging from 1 to 3. These results also incorporate the best 
practices deduced from additional analysis summarized below. The local Lax-Friedrichs flux is unstable for 
all shock positions and polynomial degrees examined. The stability of the Roe flux 18 is much improved but 
it is still unstable for large regions of shock position. Oddly, p=2 is more stable than p=l or 3; having only 
a weak instability when the shock is near the center of the element. With p=l or 3, the flux has stronger 
instabilities when the shock is near an element boundary. The HLL 19 flux is neutrally stable for all shock 
positions when p < 3. At p = 3, the HLL flux has an unstable eigenvalue only when the shock is very close 
to the element boundary. Although none of the flux choices were absolutely stable for all shock position, the 
HLL flux is considered the most robust of the three considered because its instabilities are weaker and occur 
over a smaller region of the solution space. The HLL flux is also very simple and inexpensive to implement, 
and is used in all of the numerical experiments shown later. 

DG cannot be made stable in the vicinity of a shock simply by choosing the “right” flux. Some modifica- 
tion, such as the addition of artificial viscosity, is necessary in the vicinity of a shock wave. With the goal of 
minimizing the amount of artificial viscosity that is needed, the eigenvalue analysis is used to evaluate and 
choose among several subtle alternatives within the DG implementation. 

For well resolved flows, polynomial products in the non-linear inviscid flux can be truncated to degree 
p — 1 and still retain formal order properties (e.g. as in right-hand side of Eq. 7). In practice, however, 
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(a) local Lax-Friedrichs (b) Roe (c) HLL 


Figure 1. Eigenvalue analysis of DG in the vicinity of a shock containing element. 


simulations are generally performed on the coarsest acceptable grid to minimize computational cost. Super- 
convergence has been demonstrated 14 for viscous flows with products truncated to degree p or p + 1. As 
found in the earlier study, 13 however, any truncation of the polynomial products in the vicinity of a shock 
degrades the stability of the method. Similarly, the degree of the density-weighted projection, M , used to 
approximate the rational polynomials can be limited to p — 1 and still retain the formal order properties 
of the DG method. The analysis shows that stability is improved by increasing M to p, but that there is 
generally no benefit in taking M larger than p. 

However, an exception may occur depending on how the edge flux is evaluated. The edge flux can be 
computed by either taking the edge trace of the volume flux, or by evaluating the edge flux directly using the 
edge trace of the solution. However, the eigenvalue analysis indicates that the former approach for evaluating 
the edge flux requires that the degree of the density- weighted projection be increased considerably. Evaluating 
the edge flux from the edge trace of the solution is more robust and eliminates the edge dependency of the 
degree of the density-weighted projection, M. 

Finally, the eigenvalue analysis is applied for other flow conditions. Of particular interest is the case 
of a moving shock in a flow that is supersonic both upstream and downstream of the shock. This case is 
interesting because the choice of Riemann flux becomes irrelevant as the edge flux should always be taken 
from the upstream side of the edge. For a wide range of conditions, DG was found to be stable provided the 
downstream Mach number was not close to one. 


B. Exact solutions 

Starting from a few assumptions, it is possible to construct exact steady solutions to the DG discretization of 
the Euler equations on an element containing a shock. To illustrate the process, consider the DG discretiza- 
tion given by Eq. 6 on the unit element, —1/2 < x < 1/2, with F v = 0, p = 1 and the basis b ^ = {1, x}. Let 
U(x) denote the solution within the element. Also consider the case of flow from left to right with U = Ul 
and Mach > 1 for x < —1/2, and subsonic flow corresponding to the exact Rankine-Hugoniot conditions, 
U = Ur, for x > 1/2. The first equation from Eq. 6, k = 0, enforces flux conservation and gives simply 

F(U l ,U(- 1/2)) + F(U(1/2),U r ) = 0. 

The second equation from Eq. 6, k = 1 gives 

r 1 / 2 

/ F(U(x))dx = (l/2)F([/(l/2), U R ) + 1/2) F(U l , U (—1/2)). 

J- 1/2 

The first assumption, which can be verified later, is that the solution at x = —1/2 is supersonic. Thus, for 
a “strong” numerical flux that depends only on the upstream state whenever the solutions on both sides of 
an edge are supersonic (e.g. HLL and Roe fluxes) we have 

F(U(1/2),U r ) = —F(U l , U (—1/2)) = F(U l ) (8) 
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and 


1/2 


/ F(U(x))dx = F(U l ). (9) 

J- 1/2 

Eq. 8 determines the downstream state of the element solution solely as a function of the upstream flow 
conditions, regardless of the degree of the method, or the choice of higher-order basis functions. Note 
also that multiple solutions exist. If in addition to satisfying conservation and consistency conditions, the 
numerical flux also satisfies the Rankine-Hugoniot conditions (e.g. HLL and Roe fluxes), then two solutions 
are immediately evident: 

F(U (1/2), U R ) = F(U l ,Ur) = F{U r , Ur) or 17(1/2) = U L or U R . 

With the right state known, Eq. 9 can also be readily integrated and solved exactly for the p=l case 
(integrated here using the symbolic software Maple). Taking 17(1/2) = Ul gives the trivial result that flow 
in the element is constant at the upstream state, U(x) = Ul, and that the shock sits at the right-side edge 
of the element. However, taking 17(1/2) = Ur reveals two unsettling results: 1) the exact solution is not at 
all similar to the trial solution used in the eigenvalue analysis above and 2) at moderately low supersonic 
Mach numbers, M « 1.5, the solution on the upstream side of the element gives a negative pressure; a 
condition that causes most simulation codes to terminate. Figures 2(a-d) give the exact solution for cases 
with 17(1/2) = Ur, and upstream Mach numbers of 1.25, 1.5, 1.75, and 2.0. The conserved variables p and 
pe remain physical in each case, but pressure becomes negative above Mach 1.5. 

To gain insight into the significance of this second result, consider the linearized ID Euler equations in 
characteristic form which includes equations 

dq/dt + (it ± a)dq/dx = 0, 

where a = \fpPjp is the speed of sound. Letting q = ed^+ht) f or an y re al a gives (3 = — (u ± a)a, and it 
becomes clear that the right traveling wave grows exponentially when a becomes imaginary. This instability 
is given the label of a “model” instability because it is directly associated with the Euler equation, and may 
not be present in other model equations, such as Burger’s equation. Also, the instability is inherent to the 
model equation and initial conditions, not the method by which it is solved. As such, the instability can only 
be eliminated by changing the model equations, or by ensuring that the offending initial conditions cannot 
occur. Following the former approach, the exponential growth can be canceled by adding a diffusion term 
to the right hand side of the characteristic equation. 

dq/dt + (it ± a)dq/dx = vd 2 q/dx 2 with v > oT 1 \/ max(0, — 7P/ p). 

Assuming the smallest resolved wave can be approximated by ah/p « 7r, the viscosity required for stability 
can be estimated: 

v ~ max; 0. -7 P/p ). 

pn 

Any solution to the Euler equations containing a negative pressure would normally be considered non- 
physical. Most simulation codes fail catastrophically while attempting to compute the sound speed, a = 
VWIp, or similar quantity. However, it is possible to harden a code to tolerate isolated regions of negative 
pressure. Doing so improves the overall robustness of the method, and allows the artificial viscosity to 
be minimized. It will be shown that the solution accuracy generally improves as the artificial viscosity is 
reduced. 

The sound speed is typically used only in Riemann flux computations, certain boundary conditions, and 
in estimating the permissible time step. A strategy for each use is readily justified as follows. First, it 
is assumed that the primary dependent variables (p, pe and pit,;) are physically realistic such that F(U) 
is well defined. It is noted that a negative pressure can only occur on the upstream side of the element, 
because the solution on the downstream side of the element is locked by conservation (Eq. 8). Also note 
that as pressure decreases toward zero from a positive value, the Mach number becomes greater than one 
and increases towards infinity. Thus the flow on the upstream edge should continue to be interpreted as a 
supersonic inflow boundary as the pressure drops below zero (confirming the first assumption above). If the 
flow in the adjacent upstream element is also supersonic, then clearly the Riemann flux should simply return 
the value of the flux from the upstream element. This strategy can be implemented in the HLL flux simply 
by replacing the wave speed computation u ± a with u ± \J mai(( ), yP/p). Boundary conditions and time 
step estimations are treated in a similar manner. 
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(a) Machl.25 (b) Machl.5 




(c) Machl.75 (d) Mach2.0 


Figure 2. Exact solutions for DG with p=l for a shock containing element. 
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IV. Construction of Numerical Viscosity 


In Ref. 6, a spatial diffusion equation is employed to ensure the artificial viscosity has the required 
smoothness properties. However, the artificial viscosity need only be continuous to overcome the problems 
they described, and as presented below, a continuous artificial viscosity can be easily constructed by strictly 
local and algebraic means. The diffusion equation for viscosity described in the previous work plays another 
equally if not more important role; it adds robustness to the solution process by slowing down and smoothing 
changes in viscosity in response to changes in the solution. A local relaxation and bounding process provides 
the essential robustness in the algebraic construction proposed here. 

A continuous artificial viscosity is constructed for a triangular grid simply by defining a value for viscosity 
at each vertex, /i, . and fitting these with a bilinear function within each element. Construction of the vertex 
viscosity is also surprisingly straight forward. At steady state, the viscosity at any given vertex goes to 
the maximum element viscosity over all elements sharing that vertex. The element viscosity is directly 
proportional to a shock sensor, which is similar to those given in Refs. 5 and 6. However, the solution 
process involves several key mechanisms designed to provide robustness against transients in the solution 
process. The first mechanism is a bounding operator that guards against large localized overshoots in 
viscosity. The second mechanism is a temporal relaxation operator that guards against feedback that can 
lead to non-physical unsteadiness. Let the suffix v denote an arbitrary vertex, and the suffix e denote an 
element containing that vertex. The vertex viscosity is computed as follows: 


dg v /dt = uj(S v — fi v ) 
S v = max ||v e {<Se} 
S e = c he B(n, m, S e (f, k)) 


S e (f,k) = 


< f ~ fkj - fk > ' 
A fkn fk ^ 
B(n,m,g) = (1/m" + 1 / g n ) 


(10) 

( 11 ) 

(12) 

(13) 

(14) 


where u>, C and m are user prescribed coefficients, h e is a length scale, < ., . > denotes the L 2 inner product, 
and / is a test function discussed below. The operator B provides a smooth bound for the positive argument 
/, where m is the upper bound of the result, and n controls the sharpness of the bound. 

One key to obtaining a smooth numerical viscosity without a spatial diffusion equation is to ensure 
controlling parameters are smooth. The length scale, h e , is related to the mesh size and ensures the shock 
width remains generally proportional to the mesh size as the mesh is refined. However, h e does not arise out 
of some element integration or other discrete evaluation, and therefore it does not need to precisely equal 
the mesh size. In fact, the spatial variation of the length scale should be smooth even if the mesh is not. An 
abrupt change in h e would cause an equally abrupt change in the shock thickness. Because the shock is only 
marginally resolved, a change in shock thickness can have significant adverse consequences. The meshes used 
in this work are smooth by construction, and h e is obtained as \fV, where d is the dimensionality of the 
problem, and V is the element volume. However, unstructured meshes can be highly irregular containing, 
for example, small sliver elements. The element volume, or even its average, may not be sufficiently smooth 
in such cases. Instead, taking the median of adjoining element values followed by some elliptic smoothing 
may be required. Fortunately, this work need only be done once as a pre-processing step. 

The test function / is a component or function of the solution, and fk is the projection of / onto a 
lower degree polynomial degree space k < p. Previous works 5,6 used density as the test function / with 
k = p — 1. Thus f — fk contains just the highest degree orthogonal terms of density. The present work uses 
the product / = pP with k = p which has several key attributes. First, f — fk directly models and measures 
the higher order components of the non-linear terms in the inviscid flux that become large and inaccurate 
near unresolved features. It is a quadratic function of the dependent variables and can be computed exactly. 
The higher-order terms of the product depend on all of the terms in each component resulting in a smoother 
sensor. A sensor based on just the highest order terms may become small or zero at isolated points where 
the shock orientation within an element results in a high-order inflection point. Finally, the product test 
function results in a higher-order shock sensor that decays rapidly as the resolution of smooth flow features 
improves. 

Figure 3(a) shows the shock sensor as a function of shock width. These results were produced by 
evaluating S e (f,p) for a model shock profile defined by f(x) = a + &tanh[(a; — x 3 )/(h/h s )\ on a unit element 
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— 1/2 < x < 1/2, where h/h s is the ratio of element size to shock width, x s is the shock location within the 
element, and a and b are chosen to match jump conditions of a shock with a specified shock Mach number, 
M s . The value plotted is the maximum found over a sampling of shock positions. The numerical viscosity 
decays at a rate of h^ p+2 \ faster than the formal order of the method, as the resolution of the shock increases. 
Consequently, in influence of the numerical viscosity on smooth flow features also decays quickly as they are 
better resolved. Similar computations were performed on the lower-order shock sensor (/ = p and k = p — 1), 
and are shown in figure 3(b). 




(b) Density test function: / = p 


Figure 3. Shock sensor evaluated for a model shock profile with a shock Mach number of 2.0. 


The shock sensor approaches an asymptote that is independent of order p in the limit of very thin shocks. 
The asymptotic value, shown in figure 4, correlates well with the theoretical value of the shock thickness 
Reynolds number, Re; ls = p s c s h s / p s , as given by Shapiro. 21 However, this may be of little practical value 
because the actual shock profile produced in a simulation may differ considerably from the smooth model 
assumed above as the shock becomes less resolved. Estimates for the coefficient C for a given order and 
shock width can be obtained from Ch, P ~ l/(\/7 Re/i s S e ) using the values given in figure 3 for a given order 
and shock width. a The black line in figure 3 identifies values that are considered sufficiently well resolved. 
Because of the high-order nature of the shock sensor, large changes in C are required to obtain significant 
changes in the shock width. However, the coefficient could be parameterized as C = Ch, P C p+1 where C is 
an 0(1) control parameter. 

V. Numerical Test and Examples 

A. Steady Blunt Body 

Inviscid flow over the upstream half of a two-dimensional unit cylinder is chosen as a test case. All simulations 
are performed time accurately using the 3rd-order TVD-Runge-Kutta 22 method. Figures 5(a-d) show a grid 
and a series of solutions for a Mach 5 flow. The p = 0 case has no artificial viscosity added. All solutions 
with p > 0 are obtained using the same shock capturing parameters: C = 200, n=2, and m = 10 5 (effectively 
unbounded). All solutions are well behaved with no obvious overshoots or oscillations. The shock becomes 
visibly narrower as the degree is increased above p = 1. Figure 6 shows solutions on the stagnation streamline. 
The shock becomes narrower and the total enthalpy error drops considerably as the degree increases. Also, 
the magnitude of the artificial viscosity decreases as the degree is increased without any change to the shock 
capturing parameters. 

a The y/j arises because of the choice of non-dimensional variables, u r = \J P/p • \/''f is replace by 1/M S for implementations 
in which u r = Uao ■ 
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Figure 4. Dependence of product shock sensor, pP, on Mach number for p=4. The solid line is the inverse of 
the shock thickness Reynolds number, (p s c s h s / p s ) — 1 , for a perfect gas as given in Shapiro. 21 


Sensitivity of the results to the shock capturing parameter is evaluated by varying C from 6.25 to 1600 
by roughly factors of 2. Figure 7 shows solutions with p = 3 for a range of C, with and without the bounding 
operator. As expected, a lower viscosity produces narrower shocks and a lower total enthalpy error. The 
artificial viscosity is effectively self limiting. An increase in C produces an increase in viscosity, which results 
in a smoother solution. The smoother solution results in a smaller value of the shock sensor that moderates 
the increase in viscosity. Without imposing the bounding function (e.g setting to large) solutions are readily 
obtained for only a fairly narrow range of C: 50 < C < 200. Solutions for larger values of C produce large 
initial transients that can lead to failure of the flow solver. However, solutions can be obtained in some 
cases by starting from a nearby solution (next lowest degree or next lowest C), or by reducing the viscous 
relaxation factor ui. 

All failures at low values of C have been traced to the instability associated with a negative pressure. 
Applying the bounding operator improves the robustness and allows solutions to be obtained with either 
higher or lower viscosities. This mode is most useful when attempting to perform a simulation with the 
lowest possible value of viscosity. Solutions with the lowest stable artificial viscosity are found to be the 
most robust and the most accurate. These solutions may also have isolated occurrences of negative pressure, 
which when treated as described earlier, produce no measurable adverse effects. The minimum value of 
pressure is less than minus seven in the most accurate case from figure 7 (p = 3, C = 100, and to = 0.04). 
However, the pressure contour, shown in figure 8(a) shows no evidence of this. Zooming in twice, figures 
8(b) and (c), shows that the occurrences are indeed isolated, and are limited to extrema on the upstream 
edge of the element. Based on the analysis presented earlier, this type of feature is expected to have little 
effect on the remainder of the flow field. 

The last series of tests examines the grid convergence properties of the viscous shock capturing method. 
The grids used for this test are uniform in the body normal direction (no clustering in the shock region). 
Cases are identified by the number of elements around the cylinder, nx. The mesh size varies with the degree 
of the method p to allow comparison between solutions with a similar number of unknowns. The set of grids 
used for p = 3 are nx = 12, 18, 27, and 40. Figure 9 shows the coarsest grid used in this study. The HLL 
flux has been modified to be total enthalpy preserving. The p = 0 case is inviscid (no artificial viscosity); 
p = 1 and 3 both use an artificial viscosity coefficient of C = 50, and a bound, to, ranging from 0.6 to 1.5. 
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(c) P = 2 


(d) P = 3 


Figure 5. Non-dimensional pressure, P, for Mach 5 flow over cylinder, p = 0 is inviscid. All p > 0 use the same 
shock-capturing parameters: C = 200, n = 2, and m = 10 5 . 
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(a) Non-dimensional pressure (b) Normalized Enthalpy (c) Artificial viscosity 


Figure 6. Effect of varying order, p , on solutions along stagnation streamline for Mach 5 flow over cylinder. 





(a) Non-dimensional pressure 


(b) Normalized Enthalpy 


(c) Artificial viscosity 


Figure 7. Effect of varying C on solution along stagnation streamline for Mach 5 flow over cylinder, p = 3. 




(a) Non-dimensional pressure, M=5, 
p=3, C=50. 


(b) Enlargement of boxed region of fig- 
ure 8(a). 


(c) Enlargement of boxed region of fig- 
ure 8(b). 


Figure 8. Typical isolated occurrence of negative pressure for Mach 5 flow over cylinder, p = 3. 
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y 

Figure 9. Coarsest grid in mesh refinement study, nx=12. 

This bound is low enough to improve robustness during transients, but does not constrain the steady state 
viscosity. 

Figures 10(a-f) show pressure (top) and normalized total enthalpy (bottom) along the stagnation stream- 
line for the three finest grids. As seen in earlier results, the shock region becomes narrower as the degree 
of the method p increases. As expected, the shock also becomes narrower as the mesh is refined, with no 
adjustment to the shock capturing parameters required. However, p = 0 always gives the narrowest shock 
region, so there is a large penalty in switching from an inviscicl to a viscous approach that can only be 
overcome by going to very high order. The normalized total enthalpy also improves as the order is increased 
(as seen before) or as the mesh is refined, but only slowly. No effort was made to taylor the viscous terms 
of the Navier-Stokes equations to preserve total enthalpy. Figures 11 and 12 examine the total pressure 



Figure 10. Non-dimensional pressure and normalized total enthalpy on stagnation streamline for Mach 6 
cylinder, C=50, m ranges from 0.6 to 1.5. 

along the stagnation streamline and at the body. The total pressure has been renormalized to the exact 
value expected downstream of a Mach 6 normal shock. The highest order clearly gives the best solution 
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immediately downstream of the shock, figures 11(a) and (b), in spite of the large excursions in the shock 
transition region. The highest order solution is also clearly the best in the region downstream of the shock 
where the total pressure should be constant. Total pressure of the p = 1 case is slightly noisy but still fairly 
constant. Total pressure of the p = 0 case gradually rises reducing the total error; however, this behavior 
is still incorrect for this region of the flow. Convergence of the total pressure at the cylinder is shown in 
figure 12. Both p = 1 and p = 3 are converging first order; however they are both considerably more accurate 
than the p = 0 solution. 



- 1.8 - 1.6 - 1.4 - 1.2 - 1.5 - 1.4 - 1.3 - 1.2 - 1.1 - 1.4 - 1.3 - 1.2 - 1.1 

X 


(a) (b) (c) 

Figure 11. Normalized total pressure on stagnation streamline for Mach 6 cylinder, finest grids. 
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B. Turbulent Spot Formation in a Hypersonic Boundary Layer 

The next simulation models experiments recently performed in the Boeing/AFOSR Mach-6 Quiet Tunnel 
(BAM6T) at Purdue University. 23 The simulations serve to illustrate a case in which transient shock waves 
occur in an unsteady high-speed viscous flow, and the viscous shock capturing technique enables stable 
simulations to be performed with the DG method. In the experiment, illustrated in figures 13(a-c) taken 
from Ref. 23, a spark perturber introduces a disturbance which grows quickly into a non-linear turbulent 
spot. Pressure is measured on the tunnel wall at several points downstream of the spark. Typical ensemble- 
averaged results from the Purdue experiment are shown in figure 13(c). 

Little quantitative information describing the spark is available at this time. Two-dimensional axisym- 
metric simulations were performed to determine a numerical equivalent to the spark perturber. The spark 
is modeled as energy added to the energy equation, and is imposed by adding a prescribed unsteady source 
term S(x, y, t ) to the right hand side of Eq. 2. 

S(x,y,t) = Uf[l + cos(7r/2 min(l, r/r/))] 2 sin[ir min(l, t/td)] 4 , 

where aj is the amplitude of the forcing, r is the distance from the center of the source, rf defines the 
extent of the source region and td defines the duration of the forcing. Forcing parameters td = 2 • 10 -5 and 
r/ = 0.00127 are used in the simulations presented below, and a/ is varied to match the experimental results 
at x = 2.201. 

The unsteady simulations were performed starting from a shock-free steady flow solution described in 
Ref. 24. The initial residuals were subtracted from all following values to eliminate small transients that 
could arise due to differences in the discretization schemes and grid resolution. The Reynolds number is 
9.55 • 10 6 /to. The grids are nominally Cartesian grids that have been triangulated. The grid size in the 
streamwise direction is piecewise constant with a five to one refinement in the region where forcing is applied 
to the energy equation. Grid stretching is applied in the wall normal direction to place approximated half of 
the points within the boundary layer. The nominally Cartesian grid is truncated along prescribed Mach lines 
extending from the upstream and downstream boundaries. Grid resolutions are identified in the following 
discussion by stating the size of the nominal Cartesian grid before adding the embedded region and before 
discarding the regions outside the Mach lines. 

Figure 14 shows the time progression of a typical case. The transient forcing produces a shock that 
expands outward and exits the domain. The shock strength decays rapidly everywhere except at a localized 
focal point that travels outward along the Mach cone. The transient forcing also produces a disturbance that 
grows as it propagates downstream. The method runs successfully with no artificial viscosity at a low forcing 
amplitude of a/ = 5 • 10 4 ; however, the disturbance amplitude is lower than in the experiment. Pressure 
traces of the solution sampled at x = 2.201m are shown in figure 15. 

Increasing the forcing amplitude fails unless the viscous shock capturing technique is employed. Fig- 
ures 16 and 17 show results obtained with a forcing amplitude of 6 • 10 5 . Results shown in figures 17(a) 
and (b) were obtained by applying the viscous shock capturing method as described earlier. The numerical 
viscosity is largest near the shock and 10 to 100 times smaller than the physical viscosity over most of the 
domain away from the shock. However, strong gradients at the edge of the boundary layer trigger the shock 
sensor. Ideally, the shock sensor should only be triggered near the shock wave. A series of modifications 
were applied to the shock sensor to investigate the impact of the artificial viscosity that was added in regions 
away from the shock. The results shown in figures 17(c) and (d) were produced by a modified shock sensor 
that subtracts out the effect of the initial flow (denoted as “delta” in the figure legend). S e (f,k) in Eq. 13 
is replaced by: 

max[ 0, S e ((pP)(t),p) - S e ((pP)(0),p) ]. 

While this modification removes all mean flow effects, there are still noticeably high levels of artificial viscosity 
in the region near the disturbance (ss 10% of physical viscosity). A final modification, shown in figures 17(e) 
and (f), involves multiplying the shock sensor by a heuristic mask that is one along the path of the shock, 
and quickly drops to zero away from this region. Examining the pressure traces produced by these three 
cases, shown in figure 18, reveals that the solutions are all very similar in terms of phase and wavelength; 
however, the amplitude varies by about 5%. The disturbance amplitude is now similar to that seen in the 
experiments. Use of the viscous shock capturing allows for further tuning of the forcing amplitude, location, 
duration, and spatial extent as need, to more closely match the experimental results. 
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(a) Schematic diagram of the Boeing/ AFOSR Mach-6 Quiet Tunnel. 
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(b) Schematic of experimental setup with perturber and sensor locations denoted on the x-axis. 



Figure 13. Experiment configuration and typical results. 
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Figure 14. Evolution of solution produced by low amplitude transient forcing: p=3, no artificial viscosity. 
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Figure 15. Pressure traces on the tunnel wall at x=2.201m. Low amplitude transient forcing, a=5 • 10 4 ; no 
artificial viscosity. Results from 3 grids with p = 3. 



Figure 16. Pressure traces on the tunnel wall at x=2.201m. High amplitude transient forcing with viscous 
shock capturing: oif = 6 • 10 5 , C = 1, unbounded. Results from 3 grids with p = 3. 
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Figure 17. Shock viscosity with and without modifications. Pressure is shown on the left. The ratio of shock 
viscosity to physical viscosity is shown on the right. Figures (a) and (b) apply the standard viscous shock 
capturing using test function / = pP. Figures (c) and (d), shock sensor modified to remove mean flow effects. 
Figures (e) and (f) removes mean flow effects and also applies a masking function. 
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Figure 18. Pressure traces on the tunnel wall at x=2.201m showing the effect of modifying the shock viscosity. 


VI. Conclusions 

A shock capturing technique based on artificial viscosity is presented for high-order DG methods. The 
formulation of the artificial viscosity is continuous, compact, and does not require the solution of an auxiliary 
diffusion equation. The method gives robust and accurate solutions for flows with strong shocks, without 
producing post shock oscillations. Two analyses of DG applied to a shock containing element are presented. 
An eigenvalue analysis is used to evaluate alternative implementations of the DG method including: the 
choice of Riemann flux, the degree of polynomial expansion or truncation in non-linear terms, and the 
degree of projection operators. A second analysis that examines the exact solution of DG applied to the 
Euler equations reveals that negative pressures will occur resulting in an instability. The analysis also 
gives the minimum level of viscosity required to maintain stability. Numerical tests on an inviscid blunt 
body demonstrate that the solution across the shock is first order with mesh refinement, but that error is 
much lower than that of a p = 0 solution with no artificial viscosity. In comparisons between solutions 
with similar total degrees of freedom, the error in total pressure is reduced as the order is increased. The 
level of viscosity and solution error drops as the degree is increased without changing the shock capturing 
parameters. Simulations on a viscous hypersonic boundary layer demonstrate that the method performs well 
for complex unsteady flows with moving shocks that would be difficult to handle through a combination of 
order reduction and adaptive mesh refinement. 

Acknowledgments 

This work was funded by the Hypersonic and Supersonic projects of the Fundamental Aeronautics Pro- 
gram. The authors would like to thank Ms. Casper and Drs. Bereslr and Schneider for the use of their 
experimental results and illustrations in this paper. 

References 

1 Lomtev, I. and Karniadakis, G. E., “A discontinuous Galerkin method for the Navier-Stokes equations,” International 
Journal for Numerical Methods in Fluids, Vol. 29, 1999, pp. 587-603. 

2 Krivodonova, L., Xin, J., Remade, J.-F., Chevogeon, N., and Flaherty, J., “Shock detection and limiting with discontin- 
uous Galerkin methods for hyperbolic conservation laws,” Appl. Numer. Math, Vol. 48, No. 3, 2004, pp. 323-338. 

3 Hartmann, R. and Houston, P., “Adaptive Discontinuous Galerkin Finite Element Methods for the Compressible Euler 
Equations,” Journal of Computational Physics, Vol. 183, No. 2, 2002, pp. 508 - 532. 

4 Aliabadi, S., Tu, S.-Z., and Watts, M., “An Alternative to Limiter in Discontinuous Galerkin Finite Element Method For 
Simulation Of Compressible Flows,” AIAA paper 2004-0076, 42nd AIAA Aerospace Sciences Meeting and Exhibit, Reno, NV, 


21 of 22 


American Institute of Aeronautics and Astronautics 





Jan. 5-8, 2004. 

5 Persson, P.-O. and Peraire, J., “Sub-Cell Shock Capturing for Discontinuous Galerkin Methods,” AIAA paper 2006-112, 
44th AIAA Aerospace Sciences Meeting, Reno, NV, Jan. 9-12, 2006. 

6 Barter, G. E. and Darmofal, D. L., “Shock capturing with PDE-based artificial viscosity for DGFEM: Part I. Formulation,” 
Journal of Computational Physics , Vol. 229, 2010, pp. 1810-1827. 

7 Cockburn, B., Hou, S., and Shu, C.-W., “TVB Runge-Kutta Local Projection Discontinuous Galerkin Finite Element 
Method for Conservation Laws IV: The MultiDimensional Case,” Mathematics of Computation , Vol. 54, No. 190, 1990, pp. 545- 
581. 

8 Cockburn, B. and Shu, C.-W., “Runge-Kutta Discontinuous Galerkin Methods for Convection-Dominated Problems,” 
Journal of Scientific Computing , Vol. 16, No. 3, 2001, pp. 173-261. 

9 Tang, H. and Warnecke, G., “A Runge-Kutta discontinuous Galerkin method for the Euler equations,” Computers & 
Fluids , Vol. 34, 2005, pp. 375-398. 

10 Toulopoulos, I. and Ekaterinaris, J. A., “Discontinuous Galerkin Discretizations for Viscous Compressible Flows,” 5th 
GRACM International Congress on Computational Mechanics June 29 - July 1, 2005, 2005. 

11 Luo, H. and Baum, J. D., “A fast, p-Multigrid Discontinuous Galerkin Method for Compressible Flows at All Speeds,” 
AIAA paper 2006-0110, 44th AIAA Aerospace Sciences Meeting and Exhibit, Reno, NV, Jan. 9-12, 2006. 

12 Atkins, H. L. and Shu, C. W., “Quadrature-Free Implementation of Discontinuous Galerkin Method for Hyperbolic 
Equations,” AIAA Journal , Vol. 36, No. 5, 1998, pp. 775-782. 

13 Atkins, H. L., “Local Analysis of Shock Capturing Using Discontinuous Galerkin Methodology,” AIAA Paper 97-2032, 
13th AIAA Computational Fluid Dynamics Conference, Snowmass Village, CO, June 29- July 2, 1997. 

14 Atkins, H. L., “Super-convergence of Discontinuous Galerkin Method Applied to the Navier-Stokes Equations,” AIAA 
paper 2009-3787, 19th AIAA Computational Fluid Dynamics Conference, San Antonio, TX, June 22-25, 2009. 

15 Atkins, H. L. and Lockard, D. P., “A High-order Method using Unstructured Grids for Aeroacoustic Analysis of Realistic 
Aircraft Configurations,” AIAA paper 99-1945, 5th AIAA/CEAS Aeroacoustics Conference, Bellevue, WA, May 10-12, 1999. 

16 Hu, F. Q. and Atkins, H. L., “Eigensolution analysis of the discontinuous Galerkin method. Part I: One space dimensions,” 
J. Comput. Phys ., Vol. 182, 2002, pp. 516-545. 

17 Hu, F. Q. and Atkins, H. L., “Two-dimensional Wave Analysis of the Discontinuous Galerkin Method with Non-Uniform 
Grid and Boundary Conditions,” AIAA paper 2002-2514, 8th AIAA/CEAS Aeroacoustics Conference and Exhibit, Breckenridge, 
CO, June 2002. 

18 Roe, P. L., “Approximate Reimann Solver, Parameter Vectors, and Difference Schemes,” Journal of Computational 
Physics , Vol. 90, 1990, pp. 141-160. 

19 Harten, A., Lax, P. D., and van Leer, B., “On Upstream Differencing and Godunov-Type Schemes for Hyperbolic 
Conservation Laws,” SIAM Review , Vol. 25, No. 1, 1983, pp. 35-61. 

20 Toro, E., Spruce, M., and Speares, W., “Restoration of the contact surface in the HLL-Riemann solver,” Shock Waves, 
Vol. 4, 1994, pp. 25-34. 

21 Shapiro, A. H., The Dynamics and Thermodynamcis of Compressible Fluid Flow, The Ronald Press Company, 1953. 

22 Shu, C.-W., “Total-variation-diminishing time discretizations,” SIAM Journal of Scientific and Statistical Computing, 
Vol. 9, No. 6, 1988, pp. 1073-1084. 

23 Casper, K. M., Beresh, S. J., and Schneider, S. P., “Pressure Flucuations Beneath Turbulent Spots and Instability 
Wave Packets in a Hypersonic Boundary Layer,” AIAA Paper 2011-372, 49th AIAA Aerospace Sciences Meeting, Orlando FL., 
January, 2011. 

24 Greene, P. T., Eldredge, J. D., Zhong, X., and Kim, J., “A Numerical Study of Purdue’s Mach 6 Tunned with a Roughness 
Element,” AIAA paper 2009-174, 47nd AIAA Aerospace Sciences Meeting, Orlando, FL, Jan. 5-8, 2009. 


22 of 22 


American Institute of Aeronautics and Astronautics 



